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Whole brain functional connectomes hold promise for understanding human brain activity 
across a range of cognitive, developmental and pathological states. So called resting-state 
(rs) functional MRI studies have contributed to the brain being considered at a macroscopic 
scale as a set of interacting regions. Interactions are defined as correlation-based 
signal measurements driven by blood oxygenation level dependent (BOLD) contrast. 
Understanding the neurophysiological basis of these measurements is important in 
conveying useful information about brain function. Local coupling between BOLD fMRI 
and neurophysiological measurements is relatively well defined, with evidence that 
gamma (range) frequency EEG signals are the closest correlate of BOLD fMRI changes 
during cognitive processing. However, it is less clear how whole-brain network interactions 
relate during rest where lower frequency signals have been suggested to play a key 
role. Simultaneous EEG-fMRI offers the opportunity to observe brain network dynamics 
with high spatio-temporal resolution. We utilize these measurements to compare the 
connectomes derived from rs-fMRI and EEG band limited power (BLP). Merging this 
multi-modal information requires the development of an appropriate statistical framework. 
We relate the covariance matrices of the Hilbert envelope of the source localized 
EEG signal across bands to the covariance matrices derived from rs-fMRI with the 
means of statistical prediction based on sparse Canonical Correlation Analysis (sCCA). 
Subsequently, we identify the most prominent connections that contribute to this 
relationship. We compare whole-brain functional connectomes based on their geodesic 
distance to reliably estimate the performance of the prediction. The performance of 
predicting fMRI from EEG connectomes is considerably better than predicting EEG from 
fMRI across all bands, whereas the connectomes derived in low frequency EEG bands 
resemble best rs-fMRI connectivity. 



Keywords: brain connectivity, simultaneous EEG-fMRI, resting-state brain connectomes, statistical prediction, 
band limited power 



1. INTRODUCTION 

Large scale networks with correlated time courses have been 
consistently identified in the resting brain with functional 
Magnetic Resonance Imaging (fMRI) (Beckmann and Smith, 
2004; Varoquaux et al., 2010b), and electroencephalography 
(EEG) (Tagliazucchi et al., 2012) and magnetoencephalography 
(MEG) (Brookes et al., 2011a,b). Spontaneous neural fluctua- 
tions exhibit consistent correlation structures over a wide range 
of spatial and temporal scales and they constitute a prominent 
energy-consuming feature of the brain (Scholvinck et al., 2013; 
Smith et al., 2013). Several studies highlight their significance in 
modulating brain function and task efficiency (Bonnelle et al., 
2012). Furthermore, abnormalities of resting-state (rs) connec- 
tivity have been also implicated in several neurological diseases, 
including epilepsy, schizophrenia, attention deficit hyperactivity 
disorder, Alzheimers disease, stroke and traumatic brain injury 
(Zhang and Raichle, 2010). 

Multi-modal approaches and in particular combined elec- 
trophysiological measures with fMRI offer the opportunity to 



observe neurophysiological events in high temporal and spatial 
resolution. fMRI data are acquired as series of volumetric images, 
typically obtained every few seconds, that represent blood oxy- 
gen level-dependent (BOLD) contrast. This mechanism is related 
to the delivery of blood to active neuronal tissue and hence it 
allows indirect inference on brain function. This places a limit on 
the temporal resolution of neuronal fluctuations observed with 
rs-fMRI and complicates the interpretation of the estimated con- 
nectivity. On the other hand, in EEG, multiple electrodes are 
placed on the scalp to measure spontaneous electrical activity. 
Although temporal resolution of EEG is on the scale of mil- 
liseconds, the localization of the signal involves sophisticated 
algorithms and a priori models for both the source and the vol- 
ume conductor and yet it only achieves accuracy in the range of 
1-2 cm (Kaiboriboon et al., 2012). 

To fully exploit the advantages of combining multi-modal 
information, we need to understand the relationship between the 
underlying modalities as well as and their neurophysiological ori- 
gins (Laufs et al., 2008). Pioneering intracranial recordings have 
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established a link between the local BOLD signal and the under- 
ling neuronal activity (Logothetis et al., 2001; Logothetis, 2003; 
Mukamel et al, 2005; Magri et al, 2012; Chang et al, 2013). 
However, these studies do not capture the cooperative processes 
underpinning brain function that involves whole brain organi- 
zation. Furthermore, they are invasive and their application is 
limited in animals and in specific patient cohorts with neurologi- 
cal abnormalities. We are interested in examining the relationship 
of brain connectomes derived from simultaneous recordings of 
fMRI and EEG in rest. 

Specific EEG features from the scalp, such as occipital alpha 
and beta bands have been related to RSN observed with fMRI 
(Laufs et al., 2003; Moosmann et al., 2003). These studies have 
revealed networks with a large degree of commonality with rest- 
ing state networks such as the default mode and attentional 
networks. Investigating neuronal activity in different frequency 
bands has attracted considerable attention because it is hypothe- 
sized to subserve different roles and originate from anatomically 
separated but functionally related brain regions. For instance, 
band-limited gamma effects have been linked to enhanced neu- 
ral communication, while alpha oscillations have been related 
to functional inhibition (Scheeringa et al., 2011). These stud- 
ies along with studies of seed-based analysis (de Pasquale et al., 
2010; Brookes et al, 2011a) provided insight on the relationship 
of BOLD fMRI and EEG within specific networks. One major 
limitation of methodologies based on the topographic electro- 
physiological signatures of RSN is that the agreement between 
RSN observed with fMRI and EEG relies on the spatial relation- 
ship of the extracted networks (Razavi et al, 2013). This process 
depends on thresholding and it does not provide information 
about the intra-cerebral location of the EEG signal nor about the 
relationship between specific RSN connections and EEG rhythms 
(Jann et al, 2010). 

Recently, Brookes et al. derived resting state networks in a 
range of band-limited power (BLP) frequency ranges using MEG 
and investigated their relationship with the rs-fMRI (Brookes 
et al., 2011a,b). They used a beamforming source localization 
to map the MEG signal from sensor space to source/voxel 
space. Source localization provides spatial information that allows 
one to draw direct regions' correspondence across subjects. 
Subsequently, temporal independent component analysis (ICA) 
of the Hilbert envelope of the MEG signal highlighted brain net- 
works that closely resemble known rs-fMRI networks (Brookes 
et al., 2011b). This confirmed further the neurophysiological ori- 
gin of the resting-state networks that emerge in fMRI data (Smith 
etal, 2009). 

However, these comparisons were based on the spatial agree- 
ment between the temporal ICA components estimated across 
MEG frequency bands and the spatial ICA components derived 
from the analysis of rs-fMRI data (Brookes et al., 2011b). This 
approach is limited in that it uses non-simultaneous acquisition 
of MEG and fMRI data without any guarantee that differences 
in these environments (e.g., motion, auditory input) would not 
affect the outcome. Furthermore, temporal and spatial ICA can 
have diverging results, depending upon the spatiotemporal char- 
acteristics of the underlying sources. Whereas spatial agreement 
between the two maps is reassuring, further information about 



how the covariance structure between EEG and fMRI signals dif- 
fer is needed to fully understand their relationship. In particular 
knowledge of the key connections that contribute to the predic- 
tion of one connectome from the other may give insight into the 
parts of the network that are frequency specific and common to 
each modality. 

We develop a statistical framework to learn the relation- 
ship between connectomes derived from rs-fMRI and the BLP 
spectrums of simultaneous source-localized EEG recordings. To 
achieve this we relate the covariance structure of the Hilbert 
envelope of the source localized electrophysiological signal to the 
covariance matrices derived from rs-fMRI. A key methodological 
principle of this work is that the covariance structure of both the 
Hilbert envelopes of the EEG signal and the fMRI signal lie on 
a hypercone of symmetric positive matrices (SPD). In this man- 
ifold, the geodesic distance between covariance matrices can be 
estimated precisely. This provides a principled way of comparing 
multimodal weighted whole-brain networks/graphs within and 
across subjects. 

Statistical inference has been shown to be a useful tool in 
examining the relationship between brain connectivity variables 
because it establishes a link between different modalities and it 
allows the generalization of the results from a sample set to the 
general population (Deligianni et al., 2010, 2011b, 2013). We 
use statistical inference based on sparse Canonical Correlation 
Analysis (sCCA) (Witten et al, 2009a; Witten and Tibshirani, 
2009b) to link EEG and fMRI rs connectomes. Subsequently, sub- 
ject specific EEG connectomes can be predicted from previously 
unseen fMRI connectomes and vice-versa. The predicted and 
measured functional connectomes are compared based on their 
geodesic distance and a prediction error is estimated based on 
leave-one-out cross validation. This allows us to statistically assess 
the information context of fMRI and EEG brain connectomes 
across bands. 

This approach provides a rigorous multivariate statistical 
framework to quantify the importance of each connection in 
maximizing the relationship between EEG and fMRI connectivity. 
To this end, we extend the sCCA framework with the princi- 
ple of randomized Lasso (Meinshausen and Buhlmann, 2010) 
to identify the most prominent connections that contribute to 
this relationship. This assigns a probability to each connection 
to be selected, and it offers a principled way to control for false 
positives. The sCCA loadings provide a data-driven weighting 
that minimizes the influence of noisy and unrelated connections, 
which do not contribute to the relationship between EEG and 
fMRI. This also provides a quantitative assessment of the overall 
accuracy of source localization in deep-gray matter regions. 

2. MATERIALS AND METHODS 
2.1. IMAGING 

Simultaneous resting-state EEG-fMRI was acquired from 17 adult 
volunteers (11 males, 6 females, mean age: 32.84 ± 8.13 years). 
The subjects had their eyes open and were asked to remain 
awake and fixate on a white cross presented on a black back- 
ground. Subjects were asked to remain still and their head 
was immobilized using a vacuum cushion during scanning. 
Scalp EEG was recorded during the MRI scanning using a 64 
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channel MR-compatible electrode cap (BrainCap MR, Gilching, 
Germany) at native frequency of 1000 Hz. The electrodes were 
arranged according to the modified combinatorial nomenclature, 
referenced to FCz electrode. The electrocardiogram (ECG) was 
recorded, and EEG and MR scanner clocks were synchronized. 
Imaging data was acquired in a Siemens Avanto 1.5 T clinical 
scanner using a self-shielded gradient set with maximum gra- 
dient amplitude of 40mTm~ 1 and standard 12 channel head 
receiver coil. Resting-state fMRI data were acquired based on 
a T2* -weighted gradient-echo EPI sequence with 300 volumes, 
TRITE =2160/30 ms, 30 slices with thickness 3.0 mm (1mm 
gap), effective voxel size 3.3 x 3.3 x 4.0 mm, flip angle 75°, 
FOV210 x 210 x 120 mm. A Tl -weighted structural image was 
also obtained. Ethical approval has been obtained from the UCL 
Research Ethics Committee (project ID:4290/001) and informed 
consent has been obtained from all subjects. 

2.2. PREPROCESSING 

Tl-weighted images were processed with Freesurfer to obtain 
gray matter (GM) 68 cortical regions and 14 subcortical regions 
(Desikan et al., 2006) (Table SI). Comparisons between two net- 
works are easier to interpret when both are derived from the same 
set of nodes. Atlas-based parcellation allowed us to define cor- 
responding nodes in both fMRI and the source-localized EEG 
signal. We propagate the anatomical labels from Tl space to native 
fMRI space using affine registration (Modat et al., 2010) to avoid 
erroneous warping of the image due to the drop out of gradi- 
ent echo EPI images that result from local magnetic susceptibility 
effects. Anatomical labels are also propagated to MNI space, for 
the analysis of EEG, using non-rigid registration (Modat et al., 
2010). 



The first five volumes of rs-fMRI data are removed to avoid Tl 
effects and preprocessing of the functional data involves motion 
correction, high pass filtering (0.01 Hz) and spatial smoothing 
(5 mm) with FSL (Smith et al, 2004). To construct correspond- 
ing functional networks the fMRI signal is averaged across voxels 
within each GM ROI derived from the parcellation. The signal in 
WM and cerebrospinal fluid (CSF) is also averaged and along with 
the six motion parameters provided from FSL is linearly regressed 
out from the averaged time-series within each GM ROI. 

EEG was corrected offline for scanner (Allen et al., 2000) 
and cardiac pulse related artifacts (Allen et al., 1998) using 
Brain Vision Analyzer 2 (Brain Products, Gilching, Germany). 
Subsequently, it was down-sampled to 250 Hz and exported to a 
standard binary format, which is supported by SPM12b (www.fil. 
ion.ucl.ak.uk) (Friston, 2007). The pre-processed EEG signal was 
also visually reviewed and noisy channels due to low impedances 
(< 100 kOhm) were excluded from the main analysis. 

2.3. ANALYSES OF THE EEG SIGNAL 

Further analysis of the EEG signal is carried out with SPM12b. 
This involves the following steps also shown in Figure 1: 

• Bandpass filtering: The signal is filtered into five bands: S 
(l-4Hz), 9 (4-8Hz), a (8-13 Hz), ft (13-30Hz), and y 
(30-70 Hz). Phase delays are minimized by using zero-phase 
forward and reverse second order butterworth filter. Note 
that band-pass filtering is performed prior to source localiza- 
tion. Spatial resolution in beamforming is data dependent and 
thus it exhibits frequency dependent and time-variant magni- 
tude characteristics (Barnes and Hillebrand, 2003). Traditional 
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FIGURE 1 | Main steps toward deriving functional brain networks from the preprocessed EEG signal. 
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beamforming methods focus on narrow band signals because 
they approximate frequency independent of spatial selectivity. 

• Segmentation into epochs: The signal is segmented into (fMRI) 
TR epochs (2.16 s). 

• Definition of a head model: The standard template head model 
in SPM is used and the electrode positions are spatially trans- 
formed to match the template head. This provide reasonable 
co-registration of the original sensor positions to the MNI 
coordinate system of the template structural MRI image, even 
if individual subjects heads are considerably different from the 
template. 

• Definition of forward model: The three-shell boundary ele- 
ment method (BEM) model is used for forward modeling and 
the lead fields are estimated using the Sarvas formulas for each 
point on the canonical cortical mesh. 

• Source localization: EEG data is projected into source space 
using beamforming as implemented in SPM12b (Brookes et al., 
201 la, 2012). Source localization allows spatial correspondence 
across subjects and modalities. It has also the potential to 
remove signal artifacts, which cannot be explained by the scalar 
beamformer. For each GM cortical region, the EEG signal is 
projected from sensor space to points randomly drawn from 
the region, independently for each subject. The region's center 
is always included whereas the number of points is propor- 
tional to region's volume. In Figure 1, the red dots on the 
3D head model indicate the true density of random points 
drawn in cerebral cortex, which is around 0.7 points/cm 3 . 
These points have been picked randomly for each subject. Note 
that this approach of projecting the encephalography signal 
to specific brain locations has been used before to estimate 
thalamo-cortical coupling in MEG (Roux et al., 2013). 

• Estimation of Hilbert envelope: We use two approaches to esti- 
mate the EEG time series and we produce results independently 
for each case: (a) We estimate the Hilbert transform across 
the whole down-sampled time series (WTS). Therefore, con- 
nectivity matrices are estimated based on the down-sampled 
time resolution of the EEG signal, (b) The EEG time-series are 
estimated as the average of the absolute value of the Hilbert 
transform within each epoch (AWE). This results in EEG 
time-series with corresponding time samples to the fMRI time- 
series. This approach provided the best agreement with the 
fMRI signal in Brookes et al. (201 la). 

• Region average: Finally, within each region the Hilbert- 
transformed, source localized signal is averaged across the 
randomly distributed voxels to produce an EEG time-series per 
region. Note that similarly to the fMRI preprocessing, the first 
five epochs are not included in the average. 

2.4. ESTIMATION OF FUNCTIONAL BRAIN CONNECTOMES 

Once both EEG and fMRI average time-series have been esti- 
mated for each cortical region, we seek to derive the covariance 
structure of these signals. This assumes that the brain activity 
patterns are described by a Gaussian multidimensional station- 
ary process. In this case, the covariance matrix characterizes 
fully the statistical dependencies among the underlying signals 
(Sporns et al, 2000). We use the inverse covariance, normalized to 
unit diagonal to characterize functional connectivity. The inverse 



covariance, also called the precision matrix, is directly related 
to partial correlation, which provides a measure of connectivity 
strength between two regions once the influence of the others 
has been regressed out. The correlation coefficient cannot dis- 
tinguish between a direct signal transfer from node A-C from a 
signal transfer through other nodes, as for example from A to B 
to C. Partial correlation is the simplest approach in estimating 
direct connections. Furthermore, it offers a reasonable approxi- 
mation of network structure for a scale of networks of up to few 
hundred of nodes, which is what is used in practice. This simpli- 
fies the problem of associating EEG with fMRI brain connectivity. 
Hence, there is no need to consider indirect signal transfer from 
one region to another via others (Deligianni et al., 201 la). To pro- 
duce a well-conditioned, symmetric positive definite, (Synip), 
sample covariance matrix we use the shrinkage estimator (Kramer 
et al, 2009): 

A = kf+(l-k)t (1) 

where the sample covariance matrix is estimated as a convex 
linear combination of the unrestricted sample covariance matrix 
S and the estimator T, which is the identity matrix I. In this 
case, the optimal regularization parameter k e [0, 1] is deter- 
mined analytically based on the Ledoit-Wolf theorem (Ledoit and 
Wolf, 2004). This approach provides a systematic way to regular- 
ize the sample covariance matrix and it has been shown to greatly 
enhance inference of gene association networks (Schafer and 
Strimmer, 2005), where the number of variables n is much greater 
than the number of observations p. This approach allows one to 
estimate a well-conditioned covariance structure even when the 
number of connections grow quadratically with the number of 
ROIs without any prior information. 

2.5. PREDICTIVE MODEL BASED ON SPARSE CANONICAL 
CORRELATION ANALYSIS (SCCA) 

Canonical correlation analysis (CCA) is generally applied when 
one set of predictor variables X is to be related to another set 
of predicted variables Y and observations are available for both 
groups. Note that CCA is designed to deal with situations where 
the underlying variables are not statistically independent and, 
hence, they are inherently inter-correlated. The ultimate goal of 
CCA is to find two basis vectors (canonical vectors) u, v, one for 
each variable, so that the projections of X, Y onto these vectors, 
respectively are maximally linearly correlated. 

In CCA all variables from both sets are included in the fit- 
ted canonical vectors. However, for the purpose of studying 
brain connectivity, we are interested in sparse sets of associ- 
ated variables that would result in simultaneous multivariate 
dimensionality reduction and selection of the most relevant con- 
nections. Furthermore, it allows the emergence of interpretable 
links between EEG and fMRI connectivity data. Hence, we adapt 
sparse canonical correlation analysis (sCCA) to optimize the CCA 
criterion, subject to certain constrains (Witten and Tibshirani, 
2009b): 

maximize u , v « T X T Yv 
subject to :||m|| 2 < 1, ||v|| 2 < 1, ||m||i < c\, \\v\\\ < ci 
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IMIi < ci and ||v||i < c 2 represent the L\ (or lasso) penalty and 
they result in sparse canonical vectors u, v when the sparsity 
parameters C\ and c 2 , respectively, are chosen appropriately. Note 
that with u fixed, the criterion in Equation 2 is convex in v , and 
with v fixed, it is convex in u . Therefore, the objective function 
of this biconvex criterion increases in each step of an iterative 
algorithm (Witten and Tibshirani, 2009b): 

u <— argmax u M r X T Yv subject to : ||u|| 2 < 1, ||m||i < c\ 
v <— argmax v u T X T Yv subject to : ||v[| 2 < 1, ||v||i < c 2 

Here, we are interested in quantifying how well functional con- 
nectivity measured with EEG in different bands can predict fMRI 
brain connectivity and vice-versa. We use leave-one-out cross val- 
idation and thus for each subject s = 1, . . . , S, the sCCA model is 
trained based on the remaining S — 1 datasets. The number of 
components is estimated as the minimum of the ranks of the pre- 
dictor and predicted variables in CCA. The penalty values c\ , c 2 
are optimized in each cross-validation loop using an approach 
that permutes the rows of both the predictor and predicted vari- 
ables of the sCCA (Witten and Tibshirani, 2009b). Optimization 
takes place with exhaustive search on a grid of values. 

Subsequently, a subject-specific rs-fMRI connectome Y s is 
predicted from its previously unseen EEG connectome X s 
according to: 

Y s = (wX s )+Dv+ (4) 

Vice-versa a subject-specific EEG X s conenctome can be pre- 
dicted from its rs-fMRI connectome Y s based on the same sCCA 
solution of u and v vectors: 

X s = ((Y s v) T ) + D«+ (5) 

D is a diagonal matrix with the canonical correlation scores 
and + denotes the pseudoinverse. The sCCA optimization prob- 
lem being solved is symmetric in the two variables. However, 
the algorithm finds a local optimum, by first updating one, then 
updating the other criterion. Therefore, depending on the order 
of updates, the local optimum obtained might be different. We 
found that there was no practical difference when we reverse the 
optimization approach. 

Here, both X and Y are matrices with rows the vectorized 
upper or lower triangular part of the precision matrices across 
subjects. The diagonal elements of the normalized precision 
matrix are excluded since they are always ones. CCA applies to 
these elements without any further restrictions and hence there 
is no explicit guarantee the predicted precision matrix would 
be SPD. 

2.6. A METRIC TO COMPARE COVARIANCE MATRICES 

We are interested in estimating the similarity between predicted 
and estimated connectivity matrices based on a distance metric 
that quantifies differences in the space of covariance matrices. 
Precision and covariance matrices lie in the space of symmetric 
definite positive matrices T = Synip . The standard Euclidean 
distance on matrices, the Frobenius norm, does not account for 



the geometry of this space. Thus, this norm is ill-suited to quan- 
tify prediction errors. However, Syrtip can be parameterized as 
a Riemannian manifold using an intrinsic metric (Forstner and 
Moonen, 1999; Arsigny et al, 2006): 

d M (P, G) 2 = tr(logG^PG-i) 2 (6) 

This metric has been used successfully to build statistical frame- 
works of precision matrices Syrrip (Deligianni et al., 2011b). 
c2at is a distance metric, invariant to affine transformations and 
inversion, appropriate to quantify the distance between covari- 
ance matrices from biological data successfully (Mitteroecker and 
Bookstein, 2009). 

The dpj measure is applied in a leave-one-out cross-validation 
loop outside the sCCA algorithm to reliably estimate the out-of- 
sample modeling error. We have shown before that the cIai metric 
is suitable in quantifying the loss in a structured-output multi- 
variate regression predictive framework, because it accounts for 
the geometry of the output space, and it demonstrates evidence 
of statistical consistency (Deligianni et al., 2013). Since CCA is 
closely related to multivariate multiple regression analysis (Lutz 
and Eckert, 1994), we argue that cIai is appropriate to compare the 
prediction performance of different functional models of brain 
connectivity. 

2.7. IDENTIFICATION OF RELEVANT CONNECTIONS 

It is of great interest to identify which rs-fMRI connections 
are mostly related to functional connections derived in each 
EEG band. Toward this objective we concatenate the connections 
across all EEG bands to one variable X, whereas the rs-fMRI 
connectivity variable remains the same Y. We are interested in 
applying the same biconvex criterion described in Equation 3 to 
solve the sCCA problem that aims to find the parameters that 
maximize the linear relationship between X and Y, Equation 2. 
The concatenation of the connections across all EEG bands is 
advantageous because it does not require the choice of a spar- 
sity parameter for each band independently, which would hinder 
meaningful comparisons across bands. 

Subsequently, we modify the biconvex criterion in sCCA, 
Equation 3, based on the randomized Lasso principle 
(Meinshausen and Buhlmann, 2010). Therefore, Equation 3 
takes the following form: 

u <— argmax u (w x • u) T X T Yv subject to : 

IMI 2 < 1, IMIi < ci, w x 6 {1,0.5} ^ 

v <r- argmax v w T X r Y(v • wy) subject to : 
IMI 2 < 1, IMIi < c 2 , w Y e {1, 0.5} 

w x and w x are the coefficients weights chosen randomly equal 
to 0.5 or 1, as recommended by Meinshausen and Buhlmann 
(2010); Deligianni et al. (2013). The randomized sCCA criterion 
in Equation 7 is optimized several times, which is effectively a 
strategy of resampling the connectivity data. Note that c\ and c 2 
was chosen initially based on a permutation strategy and they 
remain the same through out the randomized Lasso iterations. 
The probability of selecting a connection is then given by the 
number of times the coefficient is selected over the number of 
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repetitions. This provides a principled control on thresholding 
false positives and it is a significant improvement over the stan- 
dard Lasso penalization, which does not provide any information 
on the statistical significance of the selected features. Another 
important benefit of the randomized Lasso is that it decreases the 
dependence of the selected coefficients on the initial choice of the 
sparsity parameter, c\ and cj. 

3. RESULTS 

We present results based on brain connectomes derived from the 
whole time series of EEG (WTS) as well as corresponding results 
derived based on brain connectomes estimated from averaging 
the Hilbert transformed EEG signal within each epoch (AWE). 

In Figure 2 we show the average functional connectivity matri- 
ces across subjects in fMRI and EEG S, 6, a, ft, and y bands, 
respectively. In both fMRI and EEG, the precision matrices have 
been estimated based on time-series across the whole experiment 
(WTS). Matrices are symmetric, since they reflect correlation and 
this implies that there is no directionality information. Each of 
the connectivity matrices has been estimated by averaging (mean) 
each connection across all subjects. All matrices have two dis- 
tinctive parallel lines to the diagonal that represent homologous 
inter-hemispheric connections. These are strong in both fMRI 
and EEG across all bands, whereas in EEG we also observe strong 
intra-hemispheric connections. The top row depicts the partial 
correlation within cortical regions and results in 68 x 68 matri- 
ces. The bottom row demonstrates the partial correlation within 
cortical and subcortical regions and results in 82 x 82 matrices. 
For the cortical regions, the top left matrix quadrant represents 
connections within the left hemisphere (lh), the bottom right 
represent connections within the right hemisphere (rh) and the 



remaining quadrants represent inter-hemispheric connections. In 
the bottom row, the subcortical regions have been added at the 
top left corner of the connectivity matrices. (The regions are 
given in Table SI and they are ordered similarly in their matrix 
representation.) 

In Figure 3, the connections with the 15% highest absolute 
value in Figure 2 (WTS) are shown as 3D graphs in MNI space. 
The top row shows partial correlation networks within corti- 
cal regions, whereas the bottom row shows partial correlation 
networks within cortical and subcortical regions. In rs-fMRI 
connectomes inter-hemispheric connections dominate, whereas 
across connectomes from each EEG band intra-hemispheric con- 
nections are predominant. In particular, brain regions are repre- 
sented with spheres. Their centers and radii represent the center of 
mass of each underlying region and its volume, respectively. The 
color-coding of the spheres corresponds to different parts/lobes 
of the brain. Connections above the 15% threshold are repre- 
sented as cylinders with salmon color when they are positive and 
slate-gray when they are negative. The diameter of the cylinder is 
proportional to the connection's strength, scaled independently in 
the fMRI connectome and the connectome from each EEG band. 

Figures 2, 3 demonstrate a relatively similar covariance struc- 
ture across the EEG frequency bands. To examine whether there 
is a broadband phenomenon where all frequency bands fluctuate 
together within ROIs, or whether they are minimally corre- 
lated, we plot the histograms of correlations for four subjects in 
Figure 4. Within each ROI, we estimated the correlation matrix 
(5-by-5) of the averaged time series (WTS) for each band. Here, 
we show the histograms of the off diagonal correlation elements 
across all cortical ROIs. These results show low correlation val- 
ues between bands and demonstrate that the whole-brain EEG 




FIGURE 2 | Average functional connectivity matrices across subjects in 
fMRI and EEG S, 9, a, fi, and y bands, respectively. In both fMRI and EEG, 
the precision matrix has been estimated across the whole time samples 
(WTS). (A-F) Depicts the partial correlation within cortical regions (68 x 68), 
whereas the bottom row demonstrates the partial correlation within cortical 
and subcortical regions (82 x 82). (The regions are given in Table SI and they 



are ordered similarly in their matrix representation). For the cortical regions, 
the top left matrix quadrant represents connections within left hemisphere 
(lh), the bottom right represent connections within the right hemisphere (rh), 
and the remaining quadrants represent inter-hemispheric connections. At 
(G-L), the subcortical regions have been added at the top left corner of the 
connectivity matrices. 
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FIGURE 3 | The connections with the 15% highest absolute value 
in Figure 2 (WTS) are shown as 3D graphs in standard space. 

Connections are represented as cylinders with salmon color when 
they are positive and slate-gray when they are negative. Brain regions 
are represented with spheres. Their centers and radii represent the 



centers of mass of each underlying region and its volume, 
respectively. The color-coding corresponds to different parts of the 
brain. (A-F) Shows partial correlation networks within cortical regions, 
whereas (G-L) shows partial correlation networks within cortical and 
subcortical regions. 
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FIGURE 4 | Histograms of correlations between bands across all cortical ROIs for four subjects. 



connectomes are not driven by broadband signal changes but 
rather EEG signals at different frequencies operate within the 
same networks. 

Figure 5 shows results of prediction performance and inter- 
subject variability for the case of precision matrices derived 
based on the WTS approach. Results demonstrate that sCCA 
has improved the agreement between the predicted connectivity 
matrices and the corresponding measured connectivity matri- 
ces. Note that the optimization objective of Equations 2, 3 does 
not optimize the distance between connectivity matrices directly. 
sCCA learns the relationship between EEG and fMRI connections 
across subjects and as a result the Euclidean distance between the 



predicted and measured connectomes is minimized. This usually 
results in minimizing the geodesic distance between connectomes 
too. The prediction performance are represented based on the 
cIai metric, which reflects geodesic distance between SPD matri- 
ces. The smaller the distance the more similar the connectivity 
matrices should be and subsequently the better the performance 
of the sCCA training. Figure 5A shows results based only on cor- 
tical regions that summarize the prediction performance of fMRI 
from EEG (brown box-plots), EEG from fMRI (green box-plots) 
across bands, as well as the distance between the fMRI preci- 
sion matrices and the EEG precision matrices within subjects 
(white box-plots). Figure 5B shows similar results based on both 



www.frontiersin.org 



August 2014 | Volume 8 | Article 258 | 7 



Deligianni et al 



Relating fMRI and EEG connectomes 



A 

60 



prediction of fMRI from EEG 
prediction of EEG from fMRI 
comparison of EEG and fMRI 



P 



B 



B ' - 



B B B B 



Cortical regions: Prediction performance (cIai) 




prediction of fMRI from EEG 
prediction of EEG from fMRI 
comparison of EEG and fMRI 



? H B 



B 



E 



B □ 



H " B 



□ a a b e 



6 9 a P Y 

Cortical and subcortical regions: Prediction performance (cIai) 



I a a 



B F 



Cortical regions: Inter-subject variability (d,Ai) 



fMRI 5 6 a P Y 

Cortical and subcortical regions: tnter-subject variability {cLai) 



FIGURE 5 | Results of prediction performance (WTS). This figure presents 
results of prediction performance and inter-subject variability when both fMRI 
and EEG precision matrices are estimated based on all time samples. The 
distance between the rs-fMRI precision matrices and each of the EEG 
frequency banded precision matrices estimated with cIai- The smaller the 
distance the more similar the connectivity matrices should be. (A) It shows 
results based only on cortical regions that summarize the prediction 
performance of fMRI from EEG (brown box-plots) and vice-versa (green 
box-plots) across bands, as well as the distance between the fMRI precision 



matrices and the EEG precision matrices within subjects (white box-plots), 
(B) It shows results based on both cortical and sub-cortical regions that 
summarize the prediction performance of fMRI from EEG (brown box-plots) 
and EEG from fMRI (green box-plots) across bands, as well as the distance 
between the fMRI precision matrices and the EEG precision matrices within 
subjects (white box-plots), (C) It shows inter-subject variability for the 
precision matrices estimated within cortical regions. (D) It shows 
inter-subject variability for the precision matrices estimated within cortical 
and subcortical regions. 



cortical and sub-cortical regions. In all cases, the performance of 
the predictions is estimated based on leave-one-out cross vali- 
dation. c\ and C2 have been optimized in each cross-validation 
loop according to a permutation-based algorithm [6]. The num- 
ber of components is estimated as the minimum of the ranks of 
the variables X and Y. 

The ability to predict a rs-fMRI precision matrix from an 
EEG precision matrix remains relatively similar across bands and 
it is substantially better than predicting an EEG connectivity 



matrix from a rs-fMRI precision matrix. This is also shown with 
a Wilcoxon rank-sum test, which demonstrates significant sta- 
tistical differences between the prediction performance of EEG 
from fMRI and the prediction performance of fMRI from EEG 
across all bands (p-values < le-05). On the contrary, the pre- 
diction of EEG from fMRI is considerably modulated across 
bands with the low frequency bands (S, 0, a) performing bet- 
ter, similarly to the within-subject distance between the measured 
fMRI and EEG connectomes. In Table 1 we show the p-values 
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Table 1 | P-values of Wilcoxon rank-sum test for assessing differences 
between intra-subject comparisons of EEG and fMRI across bands 
(WTS) shown in Figure 5. 

0 a p y 



(A) CORTICAL CONNECTOMES 

5 0.07 0.39 0.23 3.4e-04 

6 0.39 0.006 1.2e-05 
a 0.03 6.1e-05 
P 9.4e-03 

(B) CORTICO-SUBCORTICAL CONNECTOMES 

8 0.12 0.31 0.18 0.002 

6 0.88 0.01 0.0004 

a 0.06 0.003 

P 0.06 



Bold values indicate p < 0.05. 

of Wilcoxon rank-sum tests for assessing differences between 
intra-subject comparisons of EEG and fMRI connectomes across 
bands. Figures 5C,D shows inter-subject variability for the preci- 
sion matrices estimated within only cortical and both cortical and 
subcortical regions, respectively. Inter-subject variability in fMRI 
is considerably lower than inter-subject variability across all EEG 
bands. 

Figure 6 is similar to Figure 5 but the EEG connectivity 
matrices have been produced by averaging the Hilbert trans- 
formed signal within epochs (AWE). Therefore, each EEG time 
sample corresponds to a single fMRI time sample. (The cor- 
responding connectivity matrices and 3D graphs are shown in 
Figures SI, S2.) In this case, the distance between fMRI and 
EEG is smaller across all bands compared to Figures 5A,B. 
Nevertheless, the prediction of fMRI from EEG is better than 
the prediction of EEG from fMRI. A Wilcoxon rank-sum test 
shows significant statistical differences in 6 and bands with p- 
values of 0.04 and 0.01, respectively, for cortical connectomes and 
p-values of 0.04 and 0.005 for cortico-subcortical connectomes. 
sCCA training does not improve the performance of predicting 
EEG from fMRI compared to the original within-subject distance 
of fMRI and EEG connectivity. This may reflect the limits of the 
sCCA since d-Ai is not optimized explicitly and the original dis- 
tance of the connectivity matrices is already low. Figures 6C,D 
show inter-subject variability for only cortical and both cortical 
and subcortical regions, respectively. 

Figure 7B shows the normalized distance between the mea- 
sured EEG and fMRI precision matrices across bands when we 
use only subcortical structures, only cortical structures and both 
cortical and subcortical structures. The inclusion of the sub- 
cortical regions in the connectome increases the within subject 
distance between EEG and fMRI matrices (less similar connec- 
tomes). Several factors can account for this finding, including, 
the limitation of EEG source reconstruction in deep brain struc- 
tures. Note that d^i has been normalized based on the line fit of 
the median values of the simulation data in Figure 7A. This is 
approximately equivalent of dividing by the number of regions. 
Figure 7A shows how the cLai metric scales with the number of 
regions represented in the precision matrices. Simulation results 



come from the comparison of 1000 pairs of precision matrices 
drawn from random Whishart distributions of matrix order from 
10 to 100. To investigate further whether incorporating subcor- 
tical regions improves the prediction performance, we examined 
the performance of prediction of cortical fMRI connectomes from 
EEG cortical and cortico-subcortical connectomes. In this case, 
the number of regions in the predicted connectomes is the same 
and there is no need for any normalization. Subsequently, we 
used a paired Wilcoxon test to examine significance in each band. 
Our results showed that there is a trend that cortico-subcortical 
EEG connectomes predict cortical fMRI connectomes better than 
using cortical EEG connectomes alone. This difference is signif- 
icant in the S band (p = 0.01) and close to significance in the 6 
band (p = 0.08). 

Figure 8 demonstrate the results of 98050 randomized Lasso 
iterations for the EEG brain networks estimated based on WTS. 
(Figure S3 shows the corresponding results of the AWE case.) 
These results highlight the most prominent connections in sCCA 
from rs-fMRI (v) and EEG (u) brain connections across all bands. 
For this experiment we concatenate all the connections across all 
EEG bands to form the canonical variable X, whereas Y is the 
brain connectivity as it is measured from rs-fMRI. This allows 
us to draw the most relevant variables across all bands under 
the same sparsity parameters c\ and C2. Finally, we measure how 
many times each connection is selected out of the 98050 iterations 
and this provides us with a probability measure of confidence 
representing the importance of the underlying connection in 
maximizing the relationship between fMRI and EEG. The top row 
shows the 2% connections with the highest selection probability 
in fMRI and each EEG band for cortical regions only. The bottom 
row shows the 2% connections with the highest selection prob- 
ability for the configuration with both cortical and subcortical 
regions. We note that the selected features are mostly long-range 
connections. To our knowledge, the results of randomized Lasso 
represent the first attempt to show inter-relations between EEG 
and fMRI whole-brain connectomes. 

4. DISCUSSION 

We utilized the band-limited power envelope of the EEG signal 
to estimate an average time-series for each gray-matter cortical 
region based on a standard atlas-based parcellation. Based on 
this approach we describe functional connectivity with covari- 
ance matrices between corresponding regions across subjects and 
modalities. This allows us to compare resting-state functional 
connectivity derived across frequency bands from EEG with 
resting-state functional connectivity derived from BOLD fMRI. 
To our knowledge, we are the first to investigate the relation- 
ship between synchronous fMRI and EEG connectomes across 
frequency bands in a whole-brain, using source space analysis. 
fMRI connectivity is dominated by the inter-hemispheric con- 
nections between homologous areas, whereas brain connectivity 
derived from EEG shows a more complex pattern of connections 
composed by both intra-hemispheric and inter-hemispheric con- 
nections. We also observe that EEG connectomes in low frequency 
bands are the most similar to resting-state fMRI connectomes 
based on their geodesic distance of the underlying precision 
matrices. 
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FIGURE 6 I Results of prediction performance (AWE). This figure 
presents the same results as Figure 5 but EEG time series have been 
averaged within each epoch, which is equal to the fMRI-TR. The distance 
between the rs-fMRI precision matrices and each of the EEG frequency 
banded precision matrices estimated with cIai- The smaller the distance 
the more similar the connectivity matrices should be. (A) It shows results 
based only on cortical regions that summarize the prediction performance 
of fMRI from EEG (brown box-plots) and vice-versa (green box-plots) across 
bands, as well as the distance between the fMRI precision matrices and 



the EEG precision matrices within subjects (white box-plots), (B) It shows 
results based on both cortical and sub-cortical regions that summarize the 
prediction performance of fMRI from EEG (brown box-plots) and EEG from 
fMRI (green box-plots) across bands, as well as the distance between the 
fMRI precision matrices and the EEG precision matrices within subjects 
(white box-plots), (C) It shows inter-subject variability for the precision 
matrices estimated within cortical regions. (D) It shows inter-subject 
variability for the precision matrices estimated within cortical and 
subcortical regions. 



One possibility is that the low frequency bands in EEG 
are most predictive due to their higher signal-to-noise ratio. 
However, low frequency bands are affected from small drifts, eye 
blinks, cardiac, and respiration cycle and so on, whereas muscle 
artifacts and channels with low impedance affect higher frequen- 
cies. In addition for EEG-fMRI this is additionally complicated 
by the gradient and pulse artifacts that provide sources of struc- 
tured noise in particular in the alpha band at the slice frequency. 



Given this noise distribution, it is unlikely that the prediction dif- 
ference of the EEG bands is driven by the signal noise differences 
between bands. Nevertheless, we cannot exclude the possibility 
that differences in SNR across frequencies could explain some 
of the differences in similarity between fMRI and EEG brain 
connectomes across bands. 

Subsequently, we examine the connectivity derived from 
simultaneous EEG and fMRI by means of statistical prediction. 
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FIGURE 7 | (A) It shows how the d A! metric scales with the number of 
regions represented in the precision matrices. 1000 pairs of precision 
matrices were drawn from random Whishart distributions from matrix order 
of 10-100. Subsequently, the distance between the pair of matrices was 
estimated based on the d A i metric. (B) It shows the normalized distance 



between the measured EEG and fMRI precision matrices across bands when 
we use only subcortical structures, only cortical structures and both cortical 
and subcortical structures. d M has been normalized based on the line fit of 
the median values of the simulation data. This is approximately equivalent of 
dividing by the number of regions. 
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FIGURE 8 | Results derived from randomized Lasso for the WTS derived based on cortical regions only, whereas (G-L) shows the 

case. These reflect the 2% connections selected more often over all results obtained with precision matrices that include both cortical and 
sCCA repetitions. (A-F) Shows results with the precision matrices subcortical regions. 
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An advantage of a predictive framework of EEG and fMRI con- 
nectomes is that it removes noise that it is present in one modality 
and not the other. We use sCCA to predict EEG brain connectivity 
from fMRI and vice-versa. To evaluate the prediction perfor- 
mance we use leave-one-out cross validation and we compare 
the predicted connectivity matrix with the observed connectiv- 
ity matrix. We demonstrate that the performance of predicting 
fMRI connectivity from EEG is considerably better than pre- 
dicting EEG from fMRI across all bands. In fact, the prediction 
performance of EEG from fMRI follows a similar pattern to the 
distance between the original precision matrices, whereas the pre- 
diction performance of EEG from fMRI is relatively stable across 
bands. There is no significant improvement in prediction of fMRI 
from EEG using the joint information across multiple EEG fre- 
quency bands. Note that increasing the number of variables does 
not necessarily increase the prediction performance, since we use 
cross-validation loop to control for over fitting. 

This finding has several important implications. Firstly, it 
shows that there are signatures of rs-fMRI dynamics across EEG 
frequencies. This is consistent with the concept of nested oscil- 
lations and cross spectral coupling often found within EEG 
(Penny et al., 2008). Note that we have used envelope correla- 
tion amplitude and thus the phase information is not preserved. 
Nevertheless, if the phase-amplitude locking, which indicates 
nested oscillations, is intermittent then a large overall ampli- 
tude correlation is also expected (Penny et al, 2008). Secondly, 
it likely reflects the greater dynamic information content cap- 
tured by EEG in this particular spatial scale. Although, the spatial 
resolution of source localization is in the scale of 1-2 cm, most 
fMRI network analysis studies involve averaging the hemody- 
namic signal within larger regions. Our results indicate that in 
this spatial resolution the information carried in the EEG signal 
is richer than the averaged hemodynamic activity. In this con- 
text, the question of which EEG band represents best the fMRI 
is not important; any EEG band can provide similar connectivity 
information. This implies that scalp EEG can be used to provide 
similar information to resting state fMRI based connectomes at 
substantially reduced cost while providing much greater possibil- 
ities in dynamic information content. This might be because of 
the coarse brain parcellation, which limits spatial resolution to 
the size of the underling cortical regions. However, most current 
fMRI studies tend to examine connectivity at this scale. 

On the other hand, the inclusion of the subcortical regions 
results in more dissimilar fMRI and EEG connectomes even when 
we account for the difference in the number of regions. This may 
indicate that the highly complex cortico-subcortical interactions 
are not adequately captured with EEG alone. Cortico-subcortical 
interactions play an important role in regulating physiological 
rhythms that are associated with sleep or wakefulness, motor con- 
trol and so on. Furthermore, they have an eminent role in patho- 
logical conditions such as the propagation of epileptic activity in 
several epilepsy syndromes (Kahane and Depaulis, 2010; Moeller 
et al., 2013). Therefore, further investigation on how multi-modal 
data can improve the sensitivity in detecting these interactions 
both in space and time is crucial in discovering new treatments 
and understanding how brain networks work. Our multi-modal 
connectivity analysis demonstrate evidence that incorporating 



sub-cortical structures in EEG connectomes improves the predic- 
tion of cortical fMRI connectomes. Therefore, a cutoff in weights 
might be appropriate in some circumstances, but due to the 
sparsity constraints a weight should only be large enough to be 
influential if the corresponding edge is genuinely informative for 
the prediction. It should therefore not be necessary to explic- 
itly down-weight or ignore connections carrying little predictive 
information. We acknowledge that there is controversy in the abil- 
ity to detect subcortical sources with EEG source imaging alone 
(Muthuraman et al., 2014). However, Muthuraman et al. also 
showed that sources in deep gray matter structures are present in 
EEG data when segments with higher SNR are selected indicating 
a lower sensitivity of EEG to detect deep-gray matter sources com- 
pared to MEG data. Nevertheless, Plomp et al. linked event related 
potentials recorded with EEG with sources in the insula and sub- 
cortical areas such as the parahippocampus and the thalamus 
(Plomp et al., 2010). Furthermore, Moeller et al. demonstrated 
the ability of electrical source imaging in identifying deep sources 
in the thalamus and in revealing similar neuronal networks as 
with simultaneously acquired fMRI (Moeller et al., 2013). In any 
case, as we have discussed here, there is evidence in the literature 
and in our data to suggest that there may be some information in 
the scalp EEG which is attributable to deep sources. 

Furthermore, we showed that the connectomes derived in low 
frequency EEG bands {8, 9, and a) resemble best rs-fMRI con- 
nectomes. This conclusion results from estimating the precision 
matrices over the whole down-sampled EEG Hilbert-transformed 
time-series (WTS). When connectivity is estimated based on the 
average of the signal envelope within epochs (AWE), the geodesic 
distance between EEG and fMRI connectomes is smaller, reflect- 
ing the fact that averaging the EEG signal within epochs of equal 
duration to fMRI TR, approximates the rs-fMRI signal better. 
Effectively, this reduces the information content in the EEG in a 
way that better resembles the fluctuations observed in the BOLD 
signal. Furthermore, this temporal averaging of the EEG means 
that the difference in prediction performance between bands is 
smaller than the inter-subject variability within-band. 

In literature there is on-going controversy about which band in 
EEG mostly resembles rs-fMRI connectivity. Our results are con- 
sistent with de Pasquale et al. where a seed-based analysis was used 
to correlate the dorsal and default mode networks with sponta- 
neous MEG activity (de Pasquale et al., 2010) and showed that the 
band limited MEG signal in theta, alpha and beta bands is primar- 
ily related to BOLD fMRI connectivity. Similarly, Brookes et al. 
observed higher spatial agreement between resting-state fMRI 
and MEG in the a and (3 bands (Brookes et al., 201 lb). However, 
in de Pasquale et al.s MEG study inter hemispheric correlation 
between homologous regions was not observed in despite it being 
a typical feature of resting state fMRI studies, being observed in 
the first such study by Biswal et al. (1995). Furthermore, Cabral 
et al. found that the strength of correlation between brain regions 
peaks at the a and the lower end of the p frequency bands, both in 
MEG and in simulated connectivity based on coupled oscillators 
with parameters derived from structural networks (Cabral et al., 
2014). On the other hand, work in anaesthetized rats suggested 
that in the fMRI signal is mostly correlated to S band (epidu- 
ral) electrophysiological measures (Lu et al., 2007), whereas Magri 
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et al. highlighted a, f), and y bands as mostly related to BOLD 
fMRI spontaneous activity in anaesthetized monkeys. 

There is some consensus among studies that the best agree- 
ment between rs-fMRI and EEG signal is in the a frequency 
range. Our results also highlight low frequency bands, which 
could result from that both fMRI and EEG connectomes describe 
mostly long-distance connections due the relatively large vol- 
ume of the underlying regions. In fact, evidence suggests that 
the more distant two neural assemblies are, the longer the signal- 
conduction delay between them. This biases the maintenance of 
a phase relationship between the two signals over long cortical 
distances to low frequencies (Scholvinck et al., 2013). Also small 
time shifts in high frequencies cause proportionally large phase 
shifts, which limits correlations in high frequencies. EEG and 
fMRI provide measurements of whole-brain spontaneous activ- 
ity over a large range of spatial, temporal and spectral scales. Slow 
electrophysiological activity as it is derived from the envelope or 
power of a limited range of frequencies, also called band limited 
power (BLP), is of great interest for three reasons. First, changes 
occur over similar time scales as the BOLD signal. Secondly, it is 
related to large scale spontaneous oscillations observed between 
any pair of distant brain regions. Finally, they reflect intrinsic cou- 
pling modes that are closely related to structural connectivity and 
appear relatively constant across brain states (Engel et al., 2013; 
Woolrich et al, 2013). 

The main reason for mapping the sensors to source space, in 
combination with an atlas based analysis approach, is that it pro- 
vides a general framework that allows for an anatomical interpre- 
tation of the EEG data as well as a direct comparison with other 
networks derived from fMRI and Diffusion Weighted Imaging 
(DWI). This is important to allow the extension of our method- 
ology to pathological and atypical brains (Bellec et al, 2010). 
For example, in epilepsy, localizing accurately and specifically the 
epileptogenetic zones where seizures initiate is of tremendous 
importance for the surgical outcome. Current research shows 
that agreement between EEG and fMRI analysis in detecting 
the epileptogenetic zone correlates with good surgical outcome 
(Thornton et al, 2010). Our framework could be extended to 
shed light on how to interpret observations when there is no 
multi-modal agreement. For example, examining whether and 
how the relationship between fMRI and EEG brain networks 
differ in different brain states and pathological conditions is of 
particular interest in current clinical neuroscience studies. 

4.1. METHODOLOGICAL CONSIDERATIONS 

Sensor level connectivity analysis is biased by the effects of vol- 
ume conduction/field spread, since there are multiple sensors 
recording the signal from the same sources. This severely affects 
the estimation of connectivity and impedes interpretation of the 
results (Hillebrand et al., 2012). We have used a state of the art 
approach to estimate the sources from EEG recordings based 
on beamforming (Brookes et al, 2012). Although, the effect 
of field spread is not completely abolished, this approach pro- 
vides a reasonable solution and it is resilient to artifacts in EEG 
acquired during fMRI such as those due to switched magnetic 
fields gradients. Another option is to analyse the imaginary part 
of the coherence, which is robust to volume conductance (Engel 



et al., 2013). However, functional connectivity based on phase 
measurements have different interpretations than envelope based 
connectivity (Engel et al., 2013). It is more variable across brain 
states and less bound to structural connectivity. The framework 
provided here can be extended to study both the power enve- 
lope and the phase of the EEG signal that could provide valuable 
insights regarding the connectivity information across modalities. 

Our analysis assumes that functional connectivity can be ade- 
quately described as a stationary process. Most current connec- 
tivity studies assume stationarity to avoid the high complexity 
involved in modeling the dynamic signal information, which lim- 
its the ability to process connectomes with more than 10-20 
regions (Smith, 2012). Nevertheless, the extension of our frame- 
work, using for example sliding-window correlations, to examine 
the dynamic complexity of the underlying signals is of particular 
interest (Brookes et al., 2014). Here, we examine brain connec- 
tivity based on the precision matrix, which is the inverse of 
the covariance matrix and it reflects partial correlation. This is 
important to disentangle the influence of other regions on each 
pair-wise connection (Smith et al., 2013) and to allow direct com- 
parison between connectivity variables (Deligianni et al, 2011b, 
2013). Partial correlation not only is a reasonable approximation 
of direct connectivity among brain regions but compared to the 
usual correlation coefficient it is also more resilient to common 
underlying noise sources. 

The inversion of the covariance matrix requires a well- 
conditioned SPD matrix. This problem is also known as covari- 
ance selection, and in the context of brain connectivity it is 
challenging due to the problem's intrinsically high dimensional 
space, and to inter-subject variability (Varoquaux et al., 2010a). 
In fact, the empirical covariance matrix results in inaccurate esti- 
mation of the precision matrix from its inverse due to numerical 
instabilities and poor estimation of its eigen structure. Here, we 
use a shrinkage estimator (Kramer et al., 2009) based on the 
Ledoit and Wolf theorem (Ledoit and Wolf, 2004). This regu- 
larizes the estimate of the precision matrix by adding a diagonal 
matrix to the sample covariance before computing its inverse. 

Other approaches to regularizing the inverse covariance matrix 
based on shrinkage have been recently proposed (Friedman et al., 
2008) and they have been suggested in estimating connectivity 
from fMRI time series (Varoquaux et al., 2010a; Smith et al., 
2013). These approaches shrink the estimated values of the preci- 
sion matrix, so that very small values that are potentially noisy 
are forced to zero and the rest are better estimated. However, 
a major challenge is how to determine the shrinkage parameter 
(Hinne et al., 2014). This is particularly important when we com- 
pare connectivity across bands and modalities. One approach is to 
use cross-validation to choose the shrinkage parameter that best 
generalizes the estimated covariance within subjects (Pedregosa 
et al., 2011). This results in connectivity matrices with consid- 
erably different sparsity across bands and thus interpretation of 
the results is not straightforward. Other approaches hypothesize a 
structure based on prior information provided either from struc- 
tural data (Deligianni et al, 2013; Hinne et al., 2014) or using 
population priors (Varoquaux et al, 2010a) and they may intro- 
duce strong biases. Furthermore, their extension in populations 
with neurological diseases is not obvious. 
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It is important to note that the canonical correlation variables 
X and Y that represent EEG and fMRI connectivity, respectively, 
are not in the form of SPD matrices. They are produced by the 
concatenation of the vectorized upper triangular matrix of each 
precision matrix across subjects. The sCCA operates on these con- 
nectivity variables based on the lasso L\ penalty, which results in 
sparse vectors u and v. Although there is no explicit constraint 
to ensure that the prediction will be an SPD matrix, we do not 
encounter this problem when we predict fMRI from EEG con- 
nectomes. On the other hand, when we predict EEG from fMRI 
connectomes, non-SPD predictions appear on average three times 
for each cross-validation. This is worse with other approaches 
of estimating the precision matrix such as the graphical lasso 
(Friedman et al., 2008; Pedregosa et al, 2011). In this case, most 
of the predictions are not SPD and therefore we cannot proceed 
further and estimate the overall prediction performance reliably. 

We used gray matter regions derived from standard atlas-based 
parcellation, which is a common approach (Hillebrand et al., 
2012). The main advantage of this whole-brain parcellation is 
that it is well-defined in subject space and produces correspond- 
ing regions across subjects and modalities. It is well known that 
atlas-based segmentations have poor functional specialization 
and regions' sizes differ considerably from a few tens of voxels to 
thousands. This would produce differences in signal to noise ratio 
of the estimated time-series. Another approach is to use regions 
drawn from functional studies. Although, these regions are more 
functionally specialized, there is no universal agreement on how 
to produce a whole-brain representation and how to propagate it 
into subject space. We expect that more functionally specialized 
regions would improve the ability of the proposed approach to 
select relevant connections and subsequent interpretation of the 
results (Deligianni et al., 2013). 

Nevertheless, our analysis shows that strong inter-hemispheric 
connectivity between homologous regions is present in both EEG 
and fMRI connectomes. This is indicated by the lines parallel to 
the diagonal in the partial correlation matrices, Figure 2. These 
correlations emerge even though in EEG the voxel time-series to 
be averaged within a region are drawn randomly for each region 
and subject. Coupling between homologous sensory areas across 
hemispheres has been also revealed with envelope correlation 
in previous seed-based studies (Engel et al., 2013). This is also 
well established in resting-state fMRI analysis independently of 
how regions are defined (voxel based or function based) (Biswal 
et al., 2010). Furthermore, evidence shows that inter-hemispheric 
connectivity has critical significance for behavior, indicating an 
important interaction between homologous regions rather than 
an effect of averaging dissimilar signals (Carter et al., 2010). 
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Table SI | Freesurfer subcortical and cortical regions used in this work to 
define brain connectomes. 

Figure SI | Average functional connectivity matrices across subjects in 
fMRI and EEG S, 0, a, p, and y bands, respectively. The Hilbert envelope of 
the EEG signal was averaged within each epoch (AWE), which is equal to 
the fMRI-TR. Therefore, both fMRI and EEG time series have the same 
number of time samples. (A-F) Depicts the partial correlation within 
cortical regions (68 x 68), whereas (G-L) demonstrates the partial 
correlation within subcortical regions (82 x 82). The regions are given in 
Table S1 and they are ordered similarly to their matrix representation. For 
the cortical regions, the top left matrix quadrant represents connections 
within left hemisphere (Ih), the bottom right represent connections within 
the right hemisphere (rh), and the remaining quadrants represent 
inter-hemispheric connections. At (G-L), the subcortical regions have 
been added at the top left corner of the connectivity matrices. This figure 
is similar to Figure 2 with a noticeable decrease of the partial correlation 
across all EEG bands. This may be due to less time samples that would 
result in under-estimation of true connectivity due to the regularization. 

Figure S2 | The connections with the 15% highest absolute value in 
Figure SI (AWE) are shown as 3D graphs in standard space. Connections 
are represented as cylinders with salmon color when they are positive and 
slate-gray when they are negative. Brain regions are represented with 
spheres. Their centers and radii represent the centers of mass of each 
underlying region and its volume, respectively. The color-coding 
corresponds to different parts of the brain. (A-F) Shows partial correlation 
networks within cortical regions, whereas (G-L) shows partial correlation 
networks within cortical and subcortical regions. 

Figure S3 | Results derived from randomized Lasso for the AWE case. 

These reflect the 2% connections with the highest probability to be 
selected overall sCCA repetitions. (A-F) Shows results with the precision 
matrices derived based on cortical regions only, whereas (G-L) shows the 
results obtained with precision matrices that include both cortical and 
subcortical regions. 
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